clear all;
clc;
%%
nDiv = 200;
cells = nDiv*nDiv*nDiv;

%%
filename1 = 'energy_nDiv200_10e5_stdPhys.log';
filename2 = 'energy_nDiv200_10e5_livermorePhys.log';
filename3 = 'energy_nDiv200_10e5_penelopePhys.log';
file1 = fopen(filename1, 'r');
file2 = fopen(filename2, 'r');
file3 = fopen(filename3, 'r');
    
% Name xCoord yCoord zCoord eDeposit doseDeposit
input1 = fscanf(file1,'%u\t%u\t%u\t%e\t%e\n',[5 inf]);
input2 = fscanf(file2,'%u\t%u\t%u\t%e\t%e\n',[5 inf]);
input3 = fscanf(file3,'%u\t%u\t%u\t%e\t%e\n',[5 inf]);
    
% Bitte daran denken das File zu schliessen, sonst gibt es Probleme
fclose(file1);
fclose(file2);
fclose(file3);

matrix1 = zeros(nDiv, nDiv, nDiv);
matrix2 = zeros(nDiv, nDiv, nDiv);
matrix3 = zeros(nDiv, nDiv, nDiv);

for n=1:numel(input1(1,:))
    matrix1(input1(1,n)+1, input1(2,n)+1, input1(3,n)+1) = input1(4,n);
end
for n=1:numel(input2(1,:))
    matrix2(input2(1,n)+1, input2(2,n)+1, input2(3,n)+1) = input2(4,n);
end
for n=1:numel(input3(1,:))
    matrix3(input3(1,n)+1, input3(2,n)+1, input3(3,n)+1) = input3(4,n);
end
clear input1 input2 input3;
clear file1 file2 file3;

%% Mittelung
m = 5;
    matrix1m = zeros(nDiv/m, nDiv/m, nDiv/m);
    matrix2m = zeros(nDiv/m, nDiv/m, nDiv/m);
    matrix3m = zeros(nDiv/m, nDiv/m, nDiv/m);
    for k=1:nDiv/m
        indexk1 = k*m-m+1;
        indexk2 = k*m;
        for p=1:nDiv/m
            indexp1 = p*m-m+1;
            indexp2 = p*m;
            for q=1:nDiv/m
                indexq1 = q*m-m+1;
                indexq2 = q*m;
                matrix1m(k,p,q) = sum(sum(sum(matrix1(indexk1:indexk2, indexp1:indexp2, indexq1:indexq2))))/(m^3);
                matrix2m(k,p,q) = sum(sum(sum(matrix2(indexk1:indexk2, indexp1:indexp2, indexq1:indexq2))))/(m^3);
                matrix3m(k,p,q) = sum(sum(sum(matrix3(indexk1:indexk2, indexp1:indexp2, indexq1:indexq2))))/(m^3);
            end
        end
        k
    end
%%

% Figure Properties setzen
fig = figure('Position', [0, 0, 1280, 800]);
load cmap;
colormap(cmap);
set(gca, 'XLim', [1 nDiv]);
set(gca, 'YLim', [1 nDiv]);

%logmatrix = log10(matrix);
clim = [10e-4 10e2];
diffqclim = [0 1];

%for n=nDiv/4-5:3*nDiv/4+5
%for n=nDiv/2-5:nDiv/2+5
for n=nDiv/2:nDiv/2
    % Logarithmischer Plot
    laxes = subplot(2,3,1);
    log_plot = pcolor(laxes, squeeze(log10(matrix1(n,:,:))));
    set(log_plot, 'EdgeColor','none');
    axis square tight;
    set(gca, 'CLim', log10(clim));
    title(sprintf('Radiation Dosimetry of Patient - Standard Physics\nlog. Skala [%.0e %.0e], Total Energy: %.0f MeV', clim(1), clim(2), sum(sum(sum(matrix1)))));
    xlabel(n)
    
    laxes = subplot(2,3,2);
    log_plot = pcolor(laxes, squeeze(log10(matrix2(n,:,:))));
    set(log_plot, 'EdgeColor','none');
    axis square tight;
    set(gca, 'CLim', log10(clim));
    title(sprintf('Radiation Dosimetry of Patient - Livermore Physics\nlog. Skala [%.0e %.0e], Total Energy: %.0f MeV', clim(1), clim(2), sum(sum(sum(matrix2)))));
    xlabel(n)
    
    % Logarithmischer Plot (Matlab ausgetrickst)
    laxes = subplot(2,3,3);
    lin_plot = pcolor(laxes, squeeze(matrix3(n,:,:)));
    set(lin_plot, 'Visible','off');
    axis square;
    set(gca, 'CLim', clim);
    colorbar('Yscale', 'log');
    set(laxes, 'Visible', 'off');
    log_axes = axes('Position', get(laxes, 'Position'));
    log_plot = pcolor(log_axes, squeeze(log10(matrix3(n,:,:))));
    set(log_plot, 'EdgeColor','none');
    axis square tight;
    set(gca, 'CLim', log10(clim));
    title(sprintf('Radiation Dosimetry of Patient - Penelope Physics\nlog. Skala [%.0e %.0e], Total Energy: %.0f MeV', clim(1), clim(2), sum(sum(sum(matrix3)))));
    xlabel(n)    
    
    diffq1 = squeeze(abs(matrix1m(ceil(n/m),:,:)-matrix2m(ceil(n/m),:,:))./(matrix1m(ceil(n/m),:,:)+matrix2m(ceil(n/m),:,:)));
    diffq1(isnan(diffq1)) = 0;
    laxes = subplot(2,3,4);
    log_plot = pcolor(laxes, diffq1);
    set(log_plot, 'EdgeColor','none');
    axis square tight;
    set(gca, 'CLim', diffqclim);
    title(sprintf('(Standard-Livermore)/(Standard+Livermore)\n%dx%dx%d->1x1x1 Mittelung, lin. Skala: [%.1f %.1f]',m,m,m, diffqclim(1), diffqclim(2)));
    xlabel(ceil(n/m));
    
    diffq2 = squeeze(abs(matrix1m(ceil(n/m),:,:)-matrix3m(ceil(n/m),:,:))./(matrix1m(ceil(n/m),:,:)+matrix3m(ceil(n/m),:,:)));
    diffq2(isnan(diffq2)) = 0;
    laxes = subplot(2,3,5);
    log_plot = pcolor(laxes, diffq2);
    set(log_plot, 'EdgeColor','none');
    axis square tight;
    set(gca, 'CLim', diffqclim);
    title(sprintf('(Standard-Penelope)/(Standard+Penelope)\n%dx%dx%d->1x1x1 Mittelung, lin. Skala: [%.1f %.1f]',m,m,m, diffqclim(1), diffqclim(2)));
    xlabel(ceil(n/m));
    
    diffq3 = squeeze(abs(matrix2m(ceil(n/m),:,:)-matrix3m(ceil(n/m),:,:))./(matrix2m(ceil(n/m),:,:)+matrix3m(ceil(n/m),:,:)));
    diffq3(isnan(diffq3)) = 0;
    laxes = subplot(2,3,6);
    log_plot = pcolor(laxes, diffq3);
    colorbar;
    set(log_plot, 'EdgeColor','none');
    axis square tight;
    set(gca, 'CLim', diffqclim);
    title(sprintf('(Livermore-Penelope)/(Livermore+Penelope)\n%dx%dx%d->1x1x1 Mittelung, lin. Skala: [%.1f %.1f]',m,m,m, diffqclim(1), diffqclim(2)));
    xlabel(ceil(n/m));
    
    frame = getframe(fig);
end


